A RESCALING VELOCITY METHOD FOR DISSIPATIVE KINETIC EQUATIONS 

APPLICATIONS TO GRANULAR MEDIA 



FRANCIS FILBET AND THOMAS REY 



Abstract. We present a new numerical algorithm based on a relative energy scaling for collisional kinetic 
equations allowing to study numerically their long time behavior, without the usual problems related to 
the change of scales in velocity variables. It is based on the knowledge of the hydrodynamic limit of the 
model considered, but is able to compute solutions for either dilute or dense regimes. Several applications 
are presented for Boltzmann-like equations. This method is particularly efficient for numerical simulations 
of the granular gases equation with dissipative energy: it allows to study accurately the long time behavior 
of this equation and is very well suited for the study of clustering phenomena. 
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Tj" \ We are interested in numerical simulations of the long time behavior of collisional kinetic equations 
such as the Boltzmann equation for granular gases. To this end, we introduce a new technique which 
is based on the information provided by the hydrodynamic fields computed from a macroscopic model 
corresponding to the original kinetic equation. Then we simply rescale the kinetic equation according to 
these hydrodynamic quantities. The reason to do so is that the change of scales in velocity is a challenging 
^ \ numerical problem when one wants to deal with solutions to dissipative kinetic equations on a fixed grid. 
Indeed, most of the usual deterministic methods fails to capture the correct long time behavior because 
of concentration or spreading over the velocity space. 

Recently, F. Filbet and G. Russo proposed in a rescaling method for space homogeneous kinetic 
equations on a fixed grid. This idea is mainly based on the self-similar behavior of the solution to the 
kinetic equation. However for space non homogeneous case, the situation is much more complicated and 
this method cannot be applied since the transport operator and the boundary conditions break down this 
self-similar behavior. Here we propose a very simple idea based on a rescaling of the kinetic equation 
according to its hydrodynamic limit. We will see that this approach can actually be applied to various 
type of collisional kinetic equations such as elastic and inelastic Boltzmann equation (also known as the 
granular gases equation). 

Among popular methods for numerical simulations of collisional gases, the most widely used are mesh- 
less. On the one hand, the Molecular Dynamics algorithm is a deterministic method up to the random 
initialization of the particles which is valid in a dense regime [30J . This method has been used to describe 
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the apparition of shocks in supersonic sand and reproduces very accurately experimentation results |39j. 
On the other hand, Direct Simulation Monte Carlo (DSMC) methods are stochastic algorithms working 
either in dense or rarefied regimes. We refer for instance to [3] for a complete review of this topic. These 
methods are really efficient in term of computational cost since their complexity grows linearly with the 
number N of particles but they are rather inaccurate since the order of accuracy is about 0(l/yN). 
Another approach consists in a direct resolution of the Boltzmann operator on a phase space grid. For 
instance, deterministic and highly accurate methods based on a spectral discretization of the collisional 
operator have been proposed by F. Filbet, G. Naldi, L. Pareschi, G. Toscani and G. Russo in [371 ESI E3] 
for the space-homogeneous setting. Although being of complexity 0(N 2 ), they are spectrally accurate 
and then need very few points to be precise. Another spectral method, inspired of the direct fast Fourier 
transform approach of A.V. Bobylev and S. Rjasanow E] using Lagrange multiplier to improve the 
conservations was also used for space- homogeneous simulations by I. Gamba and S.H. Tharkabhushanam 
in [26] . 

Due to the large number of particles involved in the study of rarefied gas dynamics, we shall adopt a 
statistical physics' point of view, by the use of Boltzmann-like kinetic equations. For a given nonnegative 
initial condition /o, we will consider a particle distribution function f £ = f £ (t, x, v), for t > 0, x £ C M. d 
and v £ solution to the initial-boundary value problem 

at e 
f £ (0,x,v) = fo(x,v). 

where the collision operator Q is a Boltzmann-like operator, which preserves at least mass and momentum. 
This equation describes numerous models such as the Boltzmann equation for elastic and inelastic collisions 
or Fokker-Planck-Landau type equations. The parameter e > is the dimensionless Knudsen number, that 
is the ratio between the mean free path of particles before a collision and the length scale of observation. 
It measures the rarefaction of the gas: the gas is rarefied if e ~ 1 and dense if e <C 1. The open set VL is a 
bounded Lipschitz-continuous domain of M. d , which means that the model (II. ip has to be supplemented 
with boundary conditions described later. 

The article is organized as follows. We present in Section [2] the rescaling velocity method and give two 
possible choices of scaling functions, the coupled Filbet-Russo scaling and the new uncoupled macroscopic 
one. We choose to use the latter in the rest of the article, and apply it to three different collision operators 
in Section [3j the Boltzmann operator, the full granular gases operator and a simplified BGK-like granular 
gases operator. Subsequently, we introduce in Section 0] the spatial boundary conditions we have to 
impose on fh the so-called Maxwell boundary conditions, a convex combination of specular and diffusive 
reflections at the boundary. We present these conditions for the kinetic model (in classical and rescaled 
variables) and for the macroscopic equations. In Section we introduce the numerical methods we shall 
use for the discretization of each term of the problem: the collision operator, the conservative transport 
term and the system of conservation laws with source term. Finally, we present in Section \6\ some 
numerical results concerning the full and simplified granular gases models, in space homogeneous and 
inhomogeneous settings, with different geometries. 
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2. The Rescaling Velocity Method 

This section is devoted to the presentation of a new scaling for equation (jl.ip allowing to follow the 
change of scales in velocity. It is an extension to the space-dependent setting of the method first introduced 
in [21], using the relative kinetic energy as a scaling function. Moreover, the idea is very close to the 
classical one used by J. A. Carrillo and G. Toscani in [TJ] and J. A. Carrillo, M. Di Francesco and G. Toscani 
in [12] for nonlinear diffusion equations. It was also extended by J. A. Carrillo and J. Vazquez in [15] to 
show the apparition of chaotic behavior for a precise (constructive) nonlinearity : the self-similar profile 
of this equation "oscillates" between Gaussian (heat equation) and Zel'dovich-Kompaneets-Barenblatt 
(porous medium equation) profiles. 
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Let us drop for simplicity the e-dependence in /. For a given positive function oj, we introduce a new 
distribution g(t,x,£) by setting 



(2.1) 



f(t,x,v) 



Oj(t, X 



1 / 

g[t,x, 



OJ (t, x) ) 



where the function oj is assumed to be an accurate measure of the "support" or scale of the distribution 
/ in velocity variables. Then according to this scaling, the distribution g should naturally "follow" either 
the concentration or the spreading in velocity of the distribution /. 

Let us now derive the kinetic equation verified by the distribution g. Differentiating relation (|2.1|) with 
respect to time yields 

df 1 



Provided that v = oj£, one also has 

v • V x f 



dt oj d 



d 9 1 duj a- (C \ 



4-1 



VxS VxOjdiv^g) 

OJ 



Then, if / is solution to the distribution g given by (|2.ip is solution to the following equation: 



dg , „ 1 
-£+oji-V x g-- 
ot OJ 



doj 
~dt 



+ oj£- V x oj div^g) 



oj- 



where Q is such that Q(g,g) = Q(f, /)■ This equation can actually be written in the more convenient 
conservative form: 



(2.2) 



|| + div x (oj £ g) - div ? 



1 du . . . _ , 
oj at 



— Q(g,g)- 



In order to make this rescaling efficient, the main difficulty is now to choose an appropriate scaling 
function oj to define completely the distribution g. 



2.1. The Generalized Filbet-Russo Scaling. The first and natural idea follows from the work 
It consists in computing the function oj directly from the distribution function /, by setting 



(2.3) oj 



12 E, 



V d Pf 

where pf and Ef are respectively the local mass and kinetic energy of / 



\v\ 2 



(2.4) Pf := / f(v)dv, E f := / f(v)^-dv 

jR d Jm d 2 

We also define for the sequel the mean velocity field Uf and temperature Tf by 

(2.5) u f := 1 I [j(v)vdv, T f :=-L(2Ef - Pf \uf\ 2 

Pf JR d ^Pf 

Assuming that / is nonnegative, the quantity oj will then provide correct information on its support: if 
/ is concentrated, oj will be small, whereas it will be large for scattered distributions. This approach has 
been shown to be very accurate for the space- homogeneous setting in [23], but it is difficult to extend to 
our case because the definition of oj yields very restrictive constraints on the moments of g, namely 



/ g(t,x,0<% = - f g(t,x,0\Z\ 2 d£, V(t,i)el + x!l. 
JR d a Jm d 
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2.2. The Macroscopic Scaling. In order to relax this constraint on the distribution function g, we 
propose a different choice of the function oj. Moreover, to avoid a strong coupling between oj and g 
leading to nonlinear terms, we do not want to compute oj directly from the distribution function g itself 
but only from macroscopic quantities which are assumed to be close enough to the one computed from /. 
Therefore, a good candidate will be a solution to the system of macroscopic equations obtained from a 
closure of the kinetic model (|l.ip . But, as we want to develop a rescaling velocity method, what we really 
need is an accurate representation of the support of the initial condition, and not the good hydrodynamic 
limit of the model. 

To find this system of equation, we choose to assume that, in the asymptotic limit e — > 0, the distribution 
function / solution to converges towards a steady state A4 PiU) t, where p is the mass, u the velocity 
and T the temperature. Then, a first order approximation with respect to e would be 

f £ = M p , u ,t + 0(e), 
where the macroscopic quantities are solutions to |28j 

d t p + div x (pu) = 0, 

d t (pu) + div x (pu®u + pT) = 0, 

d t E + dry, (« (E + p)) = G{p, u, E). 
In this last system, E is the kinetic energy given by 

2E = dpT + p\u\ 2 , 

p > is the pressure and G£K the "energy" of the collision operator. A natural choice for oj is to set 

<"> » = {%■ 

The scaling function oj becomes independent of the original distribution but still follows the change of 
scales in velocity, provided that the closure has been done properly. Once more, we want to insist on the 
fact that this quantity does not have to be the exact "temperature" of the limit equation, but only an 
information on the support of the distribution /. The conditions to apply the method are then to have 
some good estimates (by above and below if possible) of the change of scales of the equation, through the 
temperature of its solutions. The system (|2.2p - (|2.6|) is now uncoupled which is better both for a numerical 
and mathematical point of view. Moreover, it will be easier to solve thanks to the hyperbolic structure 
of (I2T61) . 

Remark 2.1. In the case where no precise estimates of the dissipation term G are known, one possibility 
would be to close the moment equations with a Maxwellian distribution Ai PlU ,T> even if we can't compute 
explicitly the kinetic energy of the collision operator for the right hand side. This would yield the following 
uncoupled system 

d t p + div x (pu) = 0, 
d t {pu) + div x (pit ® u + pTI) = 0, 

d t E + div x {u{E + pT)) = — [ Q(M p , u , T ,M p ^ T ){t,x,v)\v\ 2 dv. 

One of the advantages of this uncoupled approach is that, as we will see in the following section, it general- 
izes the closure we used for both elastic and inelastic Boltzmann equation. We are currently investigating 
it for the case of kinetic models of flocking, such as the one presented in |25j (where momentum is not 
conserved as well). 

In the following section we give several examples of application of such a method. 



(2.6) 
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3. Applications of the Rescaling Velocity Method 

We start with the simple example of the classical Boltzmann equation and then present different 
applications to inelastic kinetic models for granular gases. 

3.1. Application to the Boltzmann Equation. The Boltzmann equation describes the behavior of 
a dilute gas of particles when the only interactions taken into account are binary elastic collisions. It is 
given by equation (jl.ip where f £ (t,x,v) is the time-dependent particle distribution function in the phase 
space and the Boltzmann collision operator Qg is a quadratic operator, local in (t, x) 

(3.1) QbUJ){v) =11 B{\v-v^cos6) [fj' - /*/] dadv*, 

JR d Js d ' 1 

where we used the shorthand / = f(v), /* = /(f*)> /' = f(v'), /* = The velocities of the colliding 

pairs («,«*) and (v', v'^) are related by 

, v + If — vJ . v + v* \v — vJ 

v = H a, = — a. 

2 2 ' * 2 2 

The collision kernel B is a non-negative function which by physical arguments of invariance only depends 

on \v — v*\ and cos 9 = u ■ a, where u = (v — v*)/\v — 

Boltzmann's collision operator has the fundamental properties of conserving mass, momentum and 

energy 

/ Q B (fJ)v(v)dv = 0, <p(v) = 1, v, \v\ 2 
and satisfies the well-known Boltzmann's ff-theorem (for the space homogeneous case) 



d_ 

dt 



f f log / dv = - f Q B (J, f) log(/) dv > 0. 



The functional — J f log / dv is the entropy of the solution. Boltzmann //-theorem implies that any 
equilibrium distribution function, i.e. any function which is a maximum of the entropy, has the form of 
a locally Maxwellian distribution [171 [28] 

(3-2) M ^ = ^^^[--^T L )^ 

where the density, velocity and temperature of the gas p, u and T are computed from the distribution 
function as in (I2.4|) -( |23]l . For further details on the physical background and derivation of the Boltzmann 
equation we refer to (J7J 07] . 

Assuming that the collision kernel is given by the variable hard spheres model, that is, 

B(\z\, cos 0) = \z\ x b(cos8) 

for A S [0, 1] and b bounded, then for a distribution g given by (|2.ip . we have using (13. ID 

Q B (f,f) = u x - d Q B (g,g) 

and the rescaled equation (|2.2p reads 



dg ,. /i (h: . . _ \ ■jJ- 



(3.3) — + div x (u>£g)- div € 



1 dov . . . _ 
uj at 



x 



e 

Moreover, applying a standard closure of the Boltzmann equation by a Maxwellian distribution (|3.2|) . we 
find that the macroscopic quantities and then cv thanks to ([2.7p . are given by 

( d t p + dw x (pu) = 0, 



(3.4) 



d t (pu) + div x (pu®u + pTl) = 0, 



I d t E + div x (u(E + pT)) = 0. 
We shall now apply the rescaling method to the granular gases equation. 
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3.2. Application to the Granular Gases Equation. A granular gas is a set of particles which interact 
by energy dissipating collisions leading to inelastic collisions. This inelasticity is characterized by a collision 
mechanics where mass and momentum are conserved and kinetic energy is dissipated. Thus, the collision 
phenomenon is a non-microreversible process. The velocities of the colliding pairs (v,v*) and (v',vl) are 
related by 

/ 1 + e . , , l + e / \ 
v = v — (u ■ a) a, = f * H — (u ■ a) a, 

where u := v — v* is the relative velocity, a G S d_1 the impact direction and e G [0, 1] the dissipation 
parameter, known as restitution coefficient. It means that the energy is dissipated in the impact direction 
only and will be chosen constant. The parametrization of post-collisional velocities can also be found 
by using some properties of the model. Indeed, it is equivalent to the conservation of impulsion and 
dissipation of energy: 

v' + v' # = v + v*, 



W\ 2 + Kl 2 - M 2 - Kl 2 = -- — —\u-a\ 2 < 0. 

This microscopic mechanism allows us to describe the granular collision operator Qx- If the nonnegative 
particles distribution function / depends only on v G M. d , the collision operator Qx(f,g) can be expressed 
in the following weak form: given a smooth test function ip, 

(3.6) / Q x (f,f)ipdv ■ = \ ( \v-v*\f*f (ip' + ip'^-ip- ip*) b(cos9)dadvdv*, 

JR d ^ JR d xR d xS d - 1 

where we have used the shorthand tp' := tp(v'), tp^ := tp(vl), ip := ip(v) and V>* := ^(v*)- 

Taking ip(v) = 1, v and \v\ 2 in (|3.6p . the relations (|3.5|) yield conservation of mass and momentum at 
the kinetic level and dissipation of kinetic energy 

/ Qx(f,f)\v\ 2 dv = -D(f). 

JR d 

where D(f) is the energy dissipation functional, given by 

(3.7) D(f) = C(l-e 2 ) // / f^v-v^dvdv* > 0, 



2 



for C an explicit, nonnegative constant depending on the cross section b and the dimension. The decrease 
of the energy, together with conservation of mass and momentum imply that the equilibria of the collision 
operator are Dirac distributions p5 u (v) where the density p and momentum pu are prescribed. This is 
the major problem for the numerical study of equation (jl.ip with deterministic algorithms because of this 
concentration in velocity variables. In that case, the rescaling method (|2.1|) becomes an essential tool to 
follow the concentration in velocity of the distribution function. 

For the hard sphere kernel (|3.6p and for a rescaled distribution function g given by (|2.ip . we have 

Qi(f,f) = oj 1 - d Qx(g,g). 

The rescaled equation (|2.2|) then reads 
(3.8) || + div x (u£g) - div € 



1 duj 
OJ ot 



-Qx{g,g). 



To complete the rescaling procedure, it remains to derive a set of equations for p and E to define the 
rescaling function ui. To this end, we will assume that the restitution coefficient e is close to 1, which 
corresponds to the weak inelasticity case, in order to write the collision operator as a perturbation of 
the classical Boltzmann operator and then write the system (|2.6|) associated to granular gases. Let us 
mention that there are alternative methods for deriving hydrodynamic equations from dissipative kinetic 
models, which can be found in the review paper of M. Bisi, J. A. Carrillo and G. Spiga [6]. Here, we 
focus on the approach developed by G. Toscani [4f| I45j . which is reminiscent of a the classical works of 
Brilliantov, Goldhirsch, Jenkins, Poschel and Zanetti [10 } 127 } IHT]: we write the granular gases operator 
as a perturbation of the elastic Boltzmann operator 

(3-9) Qx(f,f) = QbUJ) + el(f,f) + 0(e 2 ), 
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where X is a quadratic, energy dissipative nonlinear friction operator given by 

Z(f,f)( v ) = Trdiv^f / \u- a\(u- a)a f* fb(cos9)dadv*\ . 

* \JR d xS d - 1 J 

Using the expansion (|3.9p of the operator with respect to e, it was shown in that we can derive 
the macroscopic system (|2.6p associated to Q%. Assuming that the inelasticity is of the same order of the 
Knudsen number: 

(3.10) 1 - e = e, 

then, up to the second order in e, the macroscopic quantities are solution to 

d t p + div x (pu) = 0, 



(3.11) 



d t (pu) + div x (p (u ® u) + pTl) = 0, 



^ dtE + div x (u(E + pT)) = —Kdp 2 T 3 / 2 , 

for an explicit and nonnegative constant [44] . The solution of this system will give the scaling function 
u) in (l2T7j) . 

Remark 3.1. Let us notice that during the evolution if hydrodynamic quantities (p,u,T) given by the 
solution to k3.11\) and the ones obtained by computing the moments of the distribution function g solution 
to \°3.&\) are completely different, it is possible to update the scaling, that is the hydrodynamic model, 
starting with new macroscopic quantities which corresponds to moments of the distribution function g. If 
this happens it only means that the velocity grid in not well fitted with the scale of the distribution function 
g and this process corresponds to adapt the grid in the velocity space. 

3.3. Application to a Simplified Granular Gases Model. The rescaling velocity method can be also 
applied to a simplified model designed to mimic the behavior of the granular gases equation with a simpler 
collision operator. Indeed, A. Astillero and A. Santos derived in j2j [3] a very simple approximation of the 
granular gases equation, containing its main features. It consists in replacing the bilinear collision operator 
Qx by a nonlinear relaxation operator, supplemented with a drag force in the velocity space, depending 
on the moment of the particle distribution function itself, in order to match the correct macroscopic 
dissipation of kinetic energy of the granular gases equation. 

The relaxation operator we will use is the so-called Ellipsoidal Statistical BGK operator (ES-BGK) pQ. 
To this aim we first define some macroscopic quantities of the particle distribution function / such as the 
opposite of the stress tensor 

Qf(t,x) = — f (v — Uf) ® (v — Uf) f(t,x,v) dv. 

Pf JR d 

Therefore the translational temperature is related to the stress tensor as 7/ = tr(Of)/d. We finally 
introduce the corrected tensor 

T f (t,x) = [(l-u)T f I + uQf] (t,x), 

which can be viewed as a linear combination of the initial stress tensor Of and of the isotropic stress 
tensor Tfl developed by a Maxwellian distribution. The parameter — oo < v < 1 is used to modify the 
value of the Prandtl number through the formula 

< Pr = — - — < +oo for v G (-00 , 1). 
1 — v 

Thus we introduce a corrected Gaussian Q[f] defined by 

0[f\ = . exp J - 



/det(27r7» 

and the corresponding collision operator is now 

(3-12) Q s (f) =r(P f ) [G[f] - f) + CpfT 1 / 2 dW v ((v-u f )f), 
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Pure specular reflection (a = 1) 




u, 



Pure diffuse reflection (a = 0) 



Figure 1 . Schematic representation of Maxwell conditions at the boundary. 



where r is a given function depending on the kinetic pressure Pf := pt Tf and C is given by 

C 



(1 " e) T- 



The coefficient e corresponds to the restitution coefficient of the granular operator Q% and will be assumed 
equal to 1 — e as in (I3.10p . This expression yields the correct granular energy dissipation (13.71) . 

This new operator replaces the inelastic collision operator by an "equivalent" elastic operator and the 
drift term acts like an external drag force which mimics the cooling effect of the particles. In particular, 
the hydrodynamic equations associated to the rescaled equation (12. 2|) with operator (|3.12p are exactly 
given by the system (|3.11|) . The operator Q$ such that Qs{g) = Qs(f) is given by 



Qs(g) 



l 



using the shorthand T g 

Og 



r g (G[g}-g) + Cup a T]l 2 div* ((£ - u g ) g) 
t {oj 2 P g ) . Equation (|2.2|) then simply reads 



dt 



+ div x (u £ g) 



{G[g\ - g) + div ? 



1 duj 



CupgTl' 2 (Z-u g ) + -—Z + t®ZV xU j)g 



where the function u is given by (|2.7p using the closed system (|3.1ip . 



4. Boundary Conditions 

In order to define completely the mathematical problem for the kinetic equation (jl.ip . one has to 
set boundary conditions on d£l. The most simple (yet physically relevant) are the so-called Maxwell 
boundary conditions [32], where a fraction 1 — a G [0, 1] of the particles is specularly (elastically) reflected 
by a wall, whereas the remaining fraction a is thermalized, leaving the wall (moving at a velocity u w ) at 
a temperature T w > (see Figure [T]). The parameter a is called accommodation coefficient, because it 
expresses the tendency of the gas to accommodate to the state of the wall. More details on the subject can 
be found in the textbook of Cercignani [IB]. The main difficulty consists to apply appropriate boundary 
conditions on the macroscopic model which are consistent with the kinetic boundary value problem. 

Let us consider a distribution function f(t,x,v) solution to (|1.1|) for (i, x,v) G M+ x Q x M. d with Q a 
bounded Lipschitz-continuous domain of M. d . We assume the boundary dQ to have a unit outer normal 
n x for all x G dQ. Then, for an incoming velocity v G M. d (namely v • n x < 0), we set 



f(t,x,v) = aKf(t,x,v) + (1 - a)Mf(t,x,v), Mx G <9Q, v ■ n x < 0, 
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with 

f Tlf(t,x,v) = f(t,x,v - 2(n x -v)n x ), 

(4-1) 

[ Mf(t,x,v) = n(t,x)Mi >Uw ,T w (v). 
In this last expression, M.\ tUw ,T w is the wall Maxwellian 

ka t \ 1 ( \ v ~ u w\ 2 \ 

M Mw {v) := ^T^exp [~^- ) 

and // insures global mass conservation: 



(4.2) / f(t,x,v)v-n x dv = -p(t,x) Mi, Uw ,t w (y)v ■ n x dv. 

Jv-n x >0 Jv-n x <0 

Thanks to these definitions, it is classical (see for example [16]) that the following theorem holds: 

Theorem 4.1. Let f be a smooth solution to the kinetic equation (jl.ip with Maxwell boundary conditions 
m . Then, 

(1) the mass of f is globally conserved for any accommodation coefficient a £ [0,1]; 

(2) if Q preserves the kinetic energy, then in the specular case a = 1 the energy of f is globally 
conserved ; 

(3) if Q produces entropy, namely 

Qs(/,/)log(/)cfo>0 

then the entropy functional f f log / dv dx is dissipated for any accommodation coefficient a G 
[0,1]- 

We shall present the boundary conditions associated to the Maxwell conditions for the rescaled dis- 
tribution g and the hydrodynamics quantites (p,pu,E) solution to the Euler' system. We split in two 
parts, pure specular and pure diffusive boundary conditions. 

4.1. Specular Boundary Conditions. The natural translation of a specular reflection of a particle of 
velocity v in the physical space, at a point x of the boundary dVt is a specular reflection at the point x, 
but with a rescaled velocity ( = !)w. More precisely, for an incoming velocity £ (namely £ ■ n x < 0), we 
set 

Kg{t, x,Q=9 (t, x, [£ - 2(n x ■ f ) n x ]) . 
Concerning the macroscopic quantities, we need to impose at the fluid level the so-called slip boundary 
condition [43] . It corresponds to a zero-flux boundary conditions for the mean velocity: 

u ■ n x = 0. 

This condition yields in particular the global conservation of mass and kinetic energy that we observed 
at the kinetic level in Corollary 14.11 It then remains to fix a condition for the momentum flux. To mimic 
the specular reflection, we reverse the velocity at the boundary, as presented in [3]. 

4.2. Diffusive Boundary Conditions. We want to impose a velocity u w and a temperature T w at a 
point x of the boundary, using a wall maxwellian M\ )Uw ,T w in the physical space. Then, for an incoming 
velocity £ in the rescaled space, we set 

Mg(t,x,0 = p(t,x)M ~ = (0- 

In this last expression, M. ~ ^ is the "rescaled" wall Maxwellian, with u w = w u w , T w := uJ 1 T w and p 
insures the global mass conservation: 

/ g(t,x,£,)£-n x d£ = -p(t,x) / M x ~ x 

At the macroscopic level, we impose a zero-flux boundary conditions for the mean velocity: 

u ■ n x = 0. 
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This yields the global conservation of mass for the macroscopic model, corresponding to the result of 
Corollary 14.11 We also want to fix the temperature T w > of the wall. We set 

E {x) = ^p{x)T w . 

These boundary conditions allow to solve the boundary value problems (|3.4p or (|3.1ip . 

Remark 4.2. Of course, the boundary conditions we presented here for the macroscopic model are very 
approximate, and do no reflect exactly what type of conditions have to be imposed for this model to be the 
true hydrodynamic limit of the kinetic equation . An analysis of the boundary layer has to be carried 
out properly in order to do so, as was done for example by F. Coron in [TB] for the elastic Boltzmann 
equation. Nevertheless, we will see in the numerical experiments that our conditions are accurate enough 
for the rescaling velocity method to behave properly, because we don't need the exact hydrodynamic limit 
of our problem to compute the scaling function. 

5. Numerical Algorithms 

5.1. Discretization of the Problem. We will briefly present the numerical schemes we used for the 
discretization of the uncoupled models (|2.2p and (|2.6p . We need to distinguish three independent problems: 
the treatment of the free transport term, the collision term and the system of conservation laws. 

The finite volume method for a conservative transport equation. The first numerical problem 
concerns the space discretization of the conservative transport equation 



(5.1) 



-£ + div x (a(t, X) g) = 0, V (t, X) G M+ x O, 
at 



^g(0,X)=g (X), 



for a smooth function a : M+ x O — > M. d and a Lipschitz-continuous domain O C M, d . Our treatment of 
the problem is made in the framework of finite volume schemes, on Cartesian grid. This type of methods 
is particularly well suited to solve conservative transport equation. We choose to use an upwind method 
with the so-called Van Leer's flux limiter, as presented in the classical paper I Kij . 

The spectral method for the collision operator. We now focus on the numerical resolution of 

(5-2) ^L = Q(g,g), 

where Q = Qb or Q%. 

Any deterministic numerical method for the computation of the Boltzmann operator requires to work 
on a bounded velocity space. Our approach consists in adding some non physical binary collisions by 
periodizing the particle distribution function and the collision operator. This implies the loss of some 
local invariants, but a careful periodization allows at least the preservation of mass. This periodization 
is the basis of spectral methods, that we will use in our numerical simulations. 

Fourier techniques for the resolution of the Boltzmann equation have been first introduced indepen- 
dently by L. Pareschi and B. Perthame in [36] and by A. Bobylev and S. Rjasanow in [7\. It has since 
been investigated by a lot of authors. One can find numerous results about this method in the article [37] 
of L. Pareschi and G. Russo. It has been recently shown to be convergent by F. Filbet and C. Mouhot in 
[21] for the elastic case and allows spectral accuracy [llj . It has then been derived for the granular gases 
operator (13. 6p by F. Filbet, L. Pareschi and G. Toscani in [23] . 

The WENO approach for a system of conservation laws. The fluid equations (13. lip can be written 
as a system of nonlinear hyperbolic conservation laws with source term 

^ + div x F(U) = G(U), 

where U : x Q — > W 1 , F,G : R n — > W 1 and n = d + 2. There exists a large number of numerical 
methods to solve this type of problem, and we choose to use finite differences. 
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The reconstruction of the fluxes will be made by a high-order Weighted Essentially Non-Oscillatory 
( WENO) method with a 5 points stencil. Indeed, it is very well adapted to the treatment of shocks in 
fluid systems. The reader can consult the lecture notes [12] of C.W. Shu for most of the details about 
this method. 

The rescaling function ui (and its time and space derivatives) given by the macroscopic scaling (|2.7|) is 
then evolved in time thanks to the WENO fluxes T- , i . 

We shall now present numerical simulations of the rescaling velocity method for energy dissipative 
equations. 

6. Numerical Simulations 

6.1. Convergence toward a Dirac mass. In this section, we will deal with a space homogeneous 
distribution / = f(t,v), solution to the granular gases equation 



(6.1) 



df 



f(t = 0) = fo, 

with an initial condition chosen as the sum of two gaussian distributions of total temperature equal to 1: 
x 1 ( ( \v-ae\ 2 \ ( l-u + cr e| 2 \ \ w , 

Mv > : = r p {-^- ) + exp [--^ ))• v » e < 

where e = (1, . . . , 1) T and a = 1/10. With this choice, we have 
(6.2) / f (v) ( v ) dv = 

JR* \ |„|2 J 

This initial datum is far from equilibrium, for all restitution coefficients e £ [0,1]. The corresponding 
scaled distribution g = g(t,£) given by (|2.ip is then solution to the homogeneous granular equation with 
a time-dependent drift term 




(6.3) 



g{t = 0) = g . 



The initial condition go is chosen by go = fo because cj(0) = 1 according to (|6.2p . 

To study the evolution of w, we propose to compare the rescaling approach proposed by F. Filbet 
and G. Russo in [23] and the macroscopic one based on the closure of the kinetic model. In the space 
homogeneous setting, the mass p = 1 is preserved. Hence, we have 



u(t) = ^j- d E{t) 

and it is enough to compute the time evolution of the kinetic energy E to obtain uo. Then on the one hand, 
following [M] , we will study the distribution g solution to (|6.3p coupled with the exact kinetic energy 

(6-4) § = ^ I Qt(9, 9 ) 1 -^^. 

at jRd 2 

On the other hand, from the result of Subsection 13.21 we replace the ordinary differential equation (|6.4|) 
by 

(6.5) ^ = -K dP *T^, 

where dT(t) = 2E(t) because of the conservation of momentum. Let us emphasize that (|6.3p - (l6.5p is an 
uncoupled model. We will then compare it with the solution to (|6.1|) in classical variables. 
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Figure 2. Test 1 - Comparison of the first marginal distribution for coupled model 
(|6.3p - (|6.4|) (points) and uncoupled model (|6.3[) - (|6.4p (dashed line), with e = 0.8, at time 
t = 0.5, 10 and 50 for N v = 32 (left) and N v = 48 (right). 

Description of the results. Test case 1. Comparison of the coupled and uncoupled models. 

We start with a two dimensional test, designed to compare the results obtained with both coupled 
(|6.3[) - (|6.4p and uncoupled models (|6.3|) - (|6.5p . in the mildly inelastic regime e = 0.8. With this test, we 
intend to show that, even if we are not computing exactly the kinetic energy (provided that we are not in 
a quasi-elastic setting e ~ 1), the uncoupled model will yield exactly the same results than the coupled 
one even for an initial datum far from equilibrium. We choose the computational domain in velocity as 
the box T>y = [—V, V] 2 with V = 9 to avoid some possible loss of mass due to the drift term. We take 
successively N v = 32 and iV,, = 48 half Fourier modes in each direction of the velocity space. 

We present in Figure [2] the time evolution of the first marginal distribution 

S (1) (Ui) := / <7(t,ei,&K 2 . 
Jr 

We can see that both models give close results for short time in the coarse and fine grids, but some 
oscillations eventually appear in the coupled model, leading to the breakout of the scheme. This is due 
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Figure 3. Test 2 - Comparison of the first marginal of the solution to (|6.ip computed 
in classical (dashed line) and rescaled variables (solid line), with e = 0.8, at time t = 0.3, 
1 and 20. 



to the wrong computation of the kinetic energy of g when we take a sparse grid in the velocity space, 
because of the heavy tails of this distribution [33]. This phenomenon does not occur anymore when we 
refine the grid. In this case, the results given by both models are in very good agreement, especially for 
the long time behavior: the equilibrium distribution is very well captured by the new uncoupled model. 

Test case 2. Comparison of the uncoupled model with the equation solved in classical variables. 

We then compare the results of this test for the uncoupled model with a solution to equation (|6,1|) 
computed in classical variables with the spectral scheme. We take V = 8 and N v = 32 half modes in each 
directions of the classical velocity space. 

We observe in Figure that both models give the same result for short time. Then, on the one hand, 
the solution obtained by the spectral scheme in classical variables breaks down for intermediate time 
(i ~ 1) because of the concentration which yields spurious oscillations. We can see in particular the loss 
of symmetry of the solution at time t = 1. On the other hand, the rescaled variables prevent concentration 
and thus allows the spectral scheme to achieve its full accuracy: the Dirac mass is very well captured 
for large time. This is also due to the correct change of scales in the support of the scaled distribution, 
thanks to the expression of co. 

In this numerical example, we can also notice that the solution to the space homogeneous equation 
exhibits a clear self-similar behavior: it converges more quickly towards a self-similar solution made of 
two "stretched" bumps than towards the equilibrium solution of the equation, which is a Dirac mass. 
This observation is perfect agreement with the theoretical results of |34| . This is even true when we 
choose an initial condition far from equilibrium and a mild restitution coefficients, results which are up 
to our knowledge not proved yet. It is important to insist on the robustness of this rescaling technique, 
which is simply based on the temperature dissipation estimate (16. 5p , since the initial datum is far from a 
Maxwellian distribution and the restitution coefficient (e ~ 0.8) is not close to one. 

6.2. Trend to Equilibrium in the Non Homogeneous Case. We consider now the full inelastic 
Boltzmann equation in dimension d on the torus 

%+v.V x f = ^l±, x G T d , v G ~R d , 
at e 

in the physical case of hard spheres collision kernel with a constant cross section b, and successively 
periodic and specular boundary conditions in x. We first introduce the hydrodynamical fields associated 
to a kinetic distribution f(t,x,v). These are the (d + 1) fields of density p (scalar) and mean velocity 
u (vector valued) defined by the formulas (|2.4p - (l2.5p . Whenever f(t,x,v) is a smooth solution to the 
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granular gases equation with periodic boundary conditions, one has the global conservation laws for mass 
and momentum 

f(t, x, v) dx dv = 0, 



dt Jjd x 

d 



, f(t,x,v)vdxdv = 0. 

dt 7T d xR d 



Therefore, without loss of generality we shall impose 



f(t,x,v)dxdv = 1, / f(t,x,v)vdxdv = 0, 

T d xR d J7 d xR d 

where we normalized the measure on the torus: m (r d ^j = 1. These conservation laws are then enough to 
uniquely determine the stationary state of the inelastic Boltzmann equation: the normalized global Dirac 
distribution in velocity and constant in space 

(6.6) foo(v)=S (v). 

Our goal here is to investigate numerically the long time behavior of the solution /. If / is any reasonable 
solution of the inelastic Boltzmann equation, satisfying certain a priori bounds of compactness, then it is 
expected that / does converge to the global distribution function /oo as t goes to +oo. The point is that 
on the opposite to the elastic Boltzmann equation, the solution does not satisfy any entropy theorem. 
Therefore, to measure the convergence to equilibrium we compute first the evolution of the global kinetic 
energy, then the evolution of the deviation from the global mass 

C(t) = \\p(t,-)-l\\ L i. 

We present numerical results for the full granular gases equation (jl.ip , using the uncoupled rescaled 
model (|3.8p - (|3.1ip . The initial condition go(x,^) = co(0,x) d fo (x, ui(0, x) £) is set as 

|2\ 



(6.7) 



/ o 

w(0, x) ■ 



V dp (x) 

with po(x) = 1 + 0.1 cos(-7T2;) and T = [0,L]. 

Description of the results. Test case 3. Global Haff's Law. 

We start with a simple one dimensional test (in space and velocity) , intended to show that the uncoupled 
model is a good approximation of equation (jl.l|) in the space inhomogeneous case. The fact that we used 
an unidimensional velocity space is relevant, because we consider inelastic collisions. Indeed, for this 
type of collision, the dissipation of kinetic energy is enough to insure that the microscopic dynamics is 
nontrivial (which is of course not the case with elastic collisions). We take L = 1 and N x = 25 points 
in space variable, a rescaled velocity box T>y of size V = 10 and N v = 32 half Fourier modes in rescaled 
variables. We also perform computations in classical velocity variables, with a velocity box of size V = 8 
and N v = 128 half modes, in order to achieve a good accuracy and avoid oscillations in short time. The 
initial condition is given by (|6.7j) . 

We compare the evolution of the global kinetic energy E(t) of both simulations. We also plot the 
theoretical asymptotic behavior of the kinetic energy (namely, the Haff's Law |29j). given formally by 

W (1 + V2K t t) 2 

We can see in Figure H] that the Haff's Law is very well captured by the rescaling method. This is also the 
case for short time in the classical method, but due to the fast apparition of spurious oscillations (even 
with N v = 128 half modes), a strong divergence from the expected global behavior occurs in this case. 
Once more, the efficiency of the rescaling method for the long time behavior of the solution becomes clear, 
since the Haff's Law is perfectly described. 
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Figure 4. Test 3 - Comparison of the global inhomogeneous cooling process in original 
and scaled variables, with e = 0.9. 



Test case 4- Oscillations in the trend to equilibrium. 

Let us investigate more precisely the long time behavior of the numerical solution obtained from the 
rescaled method. In Figures [5] and El we present the time evolution in log-log scale of 

£{t) := [|p(t,.)-l|| L i 

for an initial condition given by (16. 7p . This quantity describes the relaxation towards global equilibrium 
of the local mass pit, x) in L 1 norm. 

Numerical investigations on the elastic Boltzmann equation showed some evidences of oscillations during 
the relaxation process, for example in the local relative entropy [22], as it was conjectured by L. Desvillettes 
and C. Villani in [19] . In the inelastic case, the situation is different. On the one hand the solution f(t) 
does not converge to a smooth steady state since foo is the Dirac mass (|6.6|) . but hopefully our rescaling 
method allows to follow the concentration in velocity. On the other hand, on the evolution of C(t) we 
still observe the apparition of oscillations, but now the oscillation frequency is constant with respect to 
the logarithmic time scale, which is a big difference with the elastic case (see Figure [5]). 

We also observe in Figure that the frequency of oscillation is different whether we consider specular 
or periodic boundary conditions. This is due to the fact that this frequency depends on the size of the 
box, and that specular conditions on a box are equivalent to periodic conditions on 2 d copies of this box, 
by symmetrization. 

Finally, we point out that the scheme breaks down when £{t) becomes too small (this is particularly 
clear for the case e = 0.05 in Figure EJ). This can be avoided by refining the space and velocity grids. 

Interpretation of the results. Here, we performed simulations on the full inelastic Boltzmann equation 
in a simplified geometry (one dimension of space and velocity, periodic and specular boundary conditions, 
fixed Knudsen number) with the spectral method to observe the evolution of the quantity C(t) and to 
check numerically if such oscillations occur. It is in fact possible to give a simple interpretation of these 
oscillations thanks to the work [20] . Since this effect is observed near the global equilibrium, which is 
homogeneous in space, one can consider that the space inhomogeneous granular gases equation is a small 
perturbation of the homogeneous equation, which exhibits a self-similar behavior and the Haff's Law. 
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Figure 5. Test 4 - Oscillations of the quantity C(t) for periodic and specular boundary conditions. 



Therefore, we rescale this equation as 

f f(t,x,v) = V(t) d h(s(t),x,V(t)v), 



(6.8) 



V{t) 



s{t) 



-- 1 + (l-e)t, 

log(l + (l-e)t) 
(1-e) 



Then, one can show after straightforward computations that the distribution h is solution to the granular 
Boltzmann equation with an anti-drift term in the usual set (see the paper [33J and the references therein) 
of self similar variables (s, x, v) := (s(t),x, V(t) v) 



(6.9) 



— + v ■ V x h + (1-e) V v {vh) = 

as e 



Due to the convergence toward an equilibrium distribution H e , we can replace the collision operator 
by its linearization C e around H e in (|6.9p . The initial condition (|6.7|) actually amounts to study thanks 
to C(t) the long time behavior of the first spatial Fourier mode of h. This behavior is formally given by 
the one of t i— >• e Xjt when t — > oo for any j 6 {0, ... ,d+ 1}, where Xj is one of the d + 2 hydrodynamic 
eigenvalues of the linearized operator C e . 

At least for weak inelasticity, we have seen in Subsection 13.21 that the inelastic collision operator can 
be understood as a perturbation of the elastic Boltzmann operator. In |4L)| . we used this expansion and 
results from |34j to perform a spectral analysis of £ e , as it has already been done for the elastic Boltzmann 
operator in [20J. This allows to compute the first terms of the Taylor expansion of the hydrodynamic 
eigenvalues Xj. We give this expansion with respect to r\ = 2ir/L and 1 — e: 



(6.10) 



A; 



irjXf +^xf 



;i-e)+0(r ? 3 ) + 0((l-e) 2 ), 



where Alj- is the first order expansion in r\ of the elastic hydrodynamic eigenvalues: 



±\li + 2> i^ ' 1 }; 
o, j € {2, d + l}. 




Figure 6. Test 4 - Time evolution of C{t) (deviation from equilibrium of the local mass) 
for different values of the restitution coefficient e. 
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Restitution coefficient e 


Damping rate B 


-(1 - e)B 


0.8 


-1.05 


0.21 


0.85 


-1.4 


0.21 


0.9 


-2.05 


0.205 


0.95 


-4.1 


0.205 



Table 1. Test 4 - Influence of the restitution coefficient e on the damping rate in the 
time evolution of the quantity C 



The second order expansion in rj is given by 



d + 2 2 



77. J'G{0,1}, 



d 



M> j'6{3,..,d + l} 1 



2(d-l) 

with the viscosity coefficient and k the heat conductivity of the Navier-Stokes limit of the elastic 
Boltzmann equation. 

The expansion (I6.10p can explain the oscillations and the damping observed in Figure El Indeed, the 
first order term in r] is imaginary and will generate oscillations during the return toward equilibrium, 
of period ^yl + 2/d. Then, the period of oscillation £ will be proportional to L^ 1 , as in [22J. In our 
case, the time has not been scaled, and then according to the definition (|6.8f) of s(t), the oscillations will 
appears in a logarithmic time. The second order term in r/ and the term in (1 — e) are nonpositive and 
then will exhibit a damping effect of magnitude rj 2 a + (1 — e)B for a, B positive, according to (|6.10p . We 
can see in Table [1] that this effect can be observed numerically: when the restitution coefficient vary the 
quantity (1 — e)B remains constant, when measured in the logarithmic scale. 

6.3. Inelastic Shear Flow Problem. The uniform shear flow problem is perhaps one of the most widely 
studied nonequilibrum problem for granular gases, both experimentally and theoretically [3]. It can be 
described as follows: consider a gas enclosed between two infinite planes located at y = — 1/2 and y = 1/2, 
in opposite relative motion with respective velocities —1/2 and 1/2 in the x-direction. The planes are not 
described as realistic thermalized, reflecting walls (i.e. Maxwell boundary conditions) but as Lee-Edwards 
boundaries: when a particle of velocity v hits one of the planes, it is reemitted through the opposite one 
with a velocity v' such that the relative velocity of the particle with respect to the plane is preserved 
before and after this exchange. In consequence, the energy of the system increases during this process. 

Thus, the particles gain energy at the boundary while other particles lose energy due to the inelasticity 
of the collisions. This problem has been solved numerically using Direct Simulation Monte Carlo methods 
in OSI]- We shall consider it using the full deterministic inelastic Boltzmann equation. The gas is taken 
initially at thermal equilibrium in the two dimensional infinite tube 

f (x,v) = ^eK P (-^pj , Vxe [-1/2,1/2] xR, \/v e R 2 . 

The evolution of the system will be induced by the inelastic collisions as well as the flow created by the 
energy input at the boundaries. 

Description of the results. Test case 5. Shear flow. 

The simulations are made thanks to the rescaled model (|3.8|) . with the granular gases collision operator 
Qx, on a simplified geometry: one dimension of space and velocity with f2 = [—1/2, 1/2]. The shear flow 
problem is mimicked by using diffusive boundary conditions: we impose constant temperature T w = 1 
and velocity u w = on the boundary x = ±1/2. The restitution coefficient of the inelastic collisions 
is e = 0.9. We take N x = 100 and 150 points for the space grid and N v = 32 half Fourier modes on a 
rescaled velocity box Vy = [— V, V] of size V = 12. We compare these results with computations made 
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Figure 7. Test 5 - Normalized density p(t, x)j p and temperature T(t, x)/T(t) of the shear 
flow problem computed in classical and rescaled variables, at time t = 0.25, 0.75 and 2. 

on the model (|1.1|> in classical variables. We take in this case N x = 100 spatial points and N v = 32 half 
Fourier modes on a velocity box T>y with V = 8. 

We present in Figure [7]the evolution of the normalized density p(t, x)/~p and temperature T(t, x)/T(t), 
where ~p and T are the space averaged quantities 

T W := 17^7 / x ) dx > P : = 77^7 / x ) dx - 
| S 2 1 Jn \ll\ Jn 

On the one hand, we observe that the computation in classical and rescaled variables agree very well for 
short time. Then, the oscillations due to the concentration in velocity variable for the classical case quickly 
break the result of this computation. On the other hand, the solution obtained in rescaled variables can 
be computed for large time, thanks to the change of scales induced by the choice (|2.7I) of uj. We also 
notice the stability of the scheme with respect of the space grid. 
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The short time behavior of this quantity is close to the one presented in [3] using DSMC simulations 
on the more realistic 2d x x 2d v geometry. Indeed, we can see that the thermalized walls heat the system 
at the boundary and the collisions cool it inside the domain. This induces very strong gradients near the 
boundaries for the macroscopic quantities. When time evolves, the temperature will eventually become 
uniform inside the domain with bigger and bigger gradients near the boundary. This large time behavior 
is different of the one described in [3] where the temperature reaches a space uniform steady state. This 
might be due to the use of a simplified geometry. Indeed, in their case, the heating of the boundary 
induced by the Lee-Edwards conditions will stop because of the uniformization of the gas. 

6.4. Cluster Formation in Granular Gases. We consider now the kinetic equation (II. ip with the 
simplified granular gases operator Qs in dimension 2 of space and velocity, on the torus T 2 = [0, l] 2 : 

dt J e y ' 

with periodic boundary conditions in x. We shall present some numerical evidences of clustering for this 
simplified granular gases model. 

The problem of cluster formation in a force-free granular media is a particularly difficult test case 
of numerical methods for granular gases simulations. It has been shown to occur numerically in both 
microscopic and macroscopic simulations. More precisely, in the classical Molecular Dynamics (MD) 
simulations of [57], clusters arise from an initially homogeneous gas of particles. It was also observed 
using MD simulations in [38] that clustering is a transient phenomenon: an initially homogeneous system 
of particles interacting by inelastic collisions will tend to form clusters and then will eventually come 
back after some time to the homogeneous cooling state. Clusters were also observed using deterministic 
numerical simulations for the macroscopic quantities, thanks to the energy dissipative Navier-Stokes model 

of pi. 

We will show that this phenomenon also occurs at the mesoscopic scale, namely using a kinetic de- 
scription. The initial condition for the problem is taken as a Maxwellian distribution of homogeneous, 
isotropic momentum and temperature and uniform density: 

f ° {X,V) = (2^7 2 " eXP ("^2~) ' V (^) gT2xR2 
with U an uniformly distributed random variable over [0, l] 2 . 

Description of the results. The simulations are made using the simplified granular gases operator Q$ 
presented in Subsection 16.41 in the rescaled model (|2.2|) - (|2.6|) . The macroscopic scaling function u, given 
by (|2.7|) . is computed thanks to the energy dissipative Euler system (|3. 1 lj) . Due to the correct energy 
dissipation, this model mimics the behavior of the kinetic granular gases equation and is then consistent 
with the MD simulations of [27] . Moreover, the ES-BGK operator used in the definition of Qs gives 
the good Prandtl number for the Navier-Stokes limit of the model, which is also consistent with the 
macroscopic model of [13] . 

Test case 6. Quasi-elastic clustering. 

We take N x = 50 points in each direction of space on the torus and N v = 32 points in each direction 
of a velocity box of size V = 7. The restitution coefficient is e = 0.99 (and then e = 0.01 thanks to the 
weak inelasticity asumption (|3.10|) ). We show in Figure [5] the density contour plot given by the kinetic 
model and by the dissipative Euler system (|3,11|) used to compute the scaling function. We see that both 
models give close results for short time, which is consistent with the fact that the hydrodynamic limit 
of the kinetic equation (jl.ip with the Fokker-Planck operator (|3. 12j) is given at the first order in e by 
Euler system (|2.6p . Then, for larger time, clusters start to appears even with a coarse grid in space. It is 
interesting to notice that the clusters remains bounded through time, even without taking into account 
the dense gas effects through a "pair correlation function" depending of a maximal packing fraction, as 
in 



Test case 7. Strongly inelastic clustering. 
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We now present a test where we take N x = 75 points in each direction of space on the torus and 
N v = 32 points in each direction of a velocity box of size V = 7. The restitution coefficient is taken as 
e = 0.75. This example corresponds to a case where the hydrodynamic closure we used does not apply to 
the model, because e is "large". Hence, we observe in Figure that the hydrodynamic density does not 
correspond at all with the kinetic one. Nevertheless, our method work well even in this case, because the 
temperature given by the macroscopic model is of the good order of magnitude compared to the one we 
obtain with the kinetic model, as shown in Figure [TU] (there is an agreement on the temperature up to 4 
digits). Hence, the scale of the distribution function in the velocity space is well resolved by the kinetic 
scheme, even if the temperature is very low (the order of magnitude at time t = 1 is about 2.5e — 2) and 
the velocity grid very coarse (only 32 points in each direction). 

6.5. Riemann Problem for the Elastic Boltzmann Equation. To finish, we consider the full kinetic 
equation with the ES-BGK operator for elastic interaction (as presented in Section T3. 31 or in pQ), in 
two dimension of space and velocity, for an infinite domain in space: 

% + v-V x f= T -M (Glfl-f) v( ^ )rfxR2 , 
at e 
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2.576^02 
2.570e-02 
2.563^02 
2.557e-02 
2.551e-02 
2.544e-02 
2.538e-02 
2.532e-02 
2.525e-02 
2.519e-02 



Figure 10. Test 7 - Contour plot of the temperature given by the simplified granular 
gases model (top) and the dissipative Euler system (bottom), at time t = 1, for e = 0.75. 



with r(p) = 0.925 pvr/2. Even if the ES-BGK operator (as well as the elastic Boltzmann operator) does 
not exhibit concentration nor spreading in the velocity variables, our rescaling velocity method can be of 
great interest if one wants to simulate some problems with a big difference of magnitude in temperature. 

Let us then consider a gas at rest in an infinite space, which is suddenly heated at one particular point. 
The change in temperature will induce a shock, which will propagates in the gas. If this change is big 
enough, the particle distribution function will be very concentrated in velocity where the gas is cold, and 
very spread where it is hot. We will then take as an initial condition the following centered Gaussian 
distribution 

fo(x,v) = - 77^7-^—7 ,, cxp 



(2^T (x)) d / 2 
where the initial temperature To verify 

T (x) 



2T (x) 



V(x,t>) £ R x 



1, if \x\ < 0.5, 
100, if \x\ > 0.5. 



The rescaled initial condition is then set as <?o("> £) = <^(0, -) d fo (■, w(0, •) £) for all £ £ with 



w(0,a?) = ^T (x). 



We presented in Figure [TT]the first marginal of this initial distribution, in both scaled and classical velocity 
variables. We can see that if one wants to solve the kinetic equation (|1.1|) in classical variables for such 
a distribution, one has to take a very large velocity space to capture the hot part of the gas, but the 
grid also has to be very thin in order to resolve the cold part, which is much more concentrated near the 
origin. Then, the number of points in the velocity space (which is two dimensional) has to be huge. 

With the rescaling velocity method, this is not a problem anymore, because as one can see in the right 
hand side part of Figure [TTl the initial condition in the rescaled space is a simple reduced Maxwellian 
distribution, uniform in velocity. All the information on the initial distribution will be contained in the 
scaling function oj. This function will be evolved thanks to the Euler system of perfect gases (|3.4p . 

Description of the results. We take N x = 75 points in each direction of space (with absorbing con- 
ditions to simulate an infinite domain) and N v = 32 points in each direction of a velocity box of size 
V = 7. We compute our numerical solution in an hydrodynamic regime, in order to compare it with the 
compressible Euler system (I3.4p . Hence, the Knudsen number is taken as e = 0.01 in order to have close 
enough macroscopic solutions. 

Test case 8. Elastic Riemann problem. 
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Figure 11. Test 8 - First marginal in space and velocity of the initial condition for a 
shock in temperature, in both classical and scaled velocity variables 

We present in Figures [T2l and [T51 the contour plot of the density, the x and y component of the velocity, 
and the temperature, at time t = 0.07 (because of the big shock, the front is evolving really fast), for 
respectively the rescaling velocity method and the Euler system. We observe that our method is able to 
compute the evolution of this problem in a sparse velocity grid, which is not the case of the solver in 
classical variables, and that our results are in very good agreement with the hydrodynamic ones. More 
particularly, the temperature profile are very close, which is one of the reasons that the rescaling method 
is working well. We just observe some slight oscillations in the kinetic case, but nothing really bad (and 
that can be overcome at a greater computational cost by using a refined grid). 
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Figure 13. Test 8 - Contour plot of the first hydrodynamic quantities for the Riemann 
problem in temperature, given by the compressible Euler system at time t = 0.7. 
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